SETUP

knitr::opts_chunk$set(echo = TRUE, cache = TRUE)

# install.packages(
#       "ape", "phyloregion", "canaper",            # phylo libraries
#       "terra", "geosphere", "prioritizr", "gdm",  # spatial libraries
#       "pals", "vegan", "tidyverse", "Rsymphony")  # other utilities

# install.packages("canaper", repos = "https://ropensci.r-universe.dev")

library(ape)
library(tidyverse)
library(terra)
library(phyloregion)
library(canaper)
library(prioritizr)
library(pals)
library(vegan)
library(geosphere)

# load additional convenience functions defined in this other script
source("../../code/R/functions.r")

select <- dplyr::select

WARMUP

Phylogenies

Let’s start with a quick primer on phylogenetic trees in R. We’ll use ape, the core R phylogenetics library. (Other libraries to be aware of for working with phylogenies include phytools, picante, and ggtree, among others.)

# load a tree and plot it
tree <- read.tree("../../data/mishler_2014/mishler_2014_acacia.tre")
plot(tree, type = "fan", show.tip.label = F)

# explore the data structure of a phylogeny
# class(tree)
# print(tree)
# str(tree)
# ?plot.phylo
# plot(tree, show.tip.label = F, type = "radial")

# plot the branches connecting a set of tips (e.g. a local community) --
# PD is the lengths of these branches (plus the edges connecting them to the root)
community <- sample(tree$tip.label, 5)
phyplot(tree, connect = community, 
        show.tip.label = F)

Exercises:

  • Load and plot the mammal phylogeny (Upham et al. 2019).
  • Plot the subtree representing the marsupials. The most recent common ancestor of the two species listed below is the crown node of all marsupials; you can find it using getMRCA() and prune it from the full tree using extract.clade().
  • Plot the full mammal tree with marsupials highlighted in red. You can use phytools::getDescendants() to index all the descendants of the marsupial MRCA.
# Load and plot the mammal phylogeny
marsupials <- c("Lestoros_inca", "Lagostrophus_fasciatus")

# Plot the subtree representing the marsupials
tree <- read.tree("../../data/upham_2019/upham_2019_tree0000.tre")
mrca <- getMRCA(tree, marsupials)
extract.clade(tree, mrca)
## 
## Phylogenetic tree with 362 tips and 361 internal nodes.
## 
## Tip labels:
##   Rhyncholestes_raphanurus, Lestoros_inca, Caenolestes_fuliginosus, Caenolestes_sangay, Caenolestes_caniventer, Caenolestes_condorensis, ...
## 
## Rooted; includes branch lengths.
# Plot the full mammal tree with marsupials highlighted in red
desc <- phytools::getDescendants(tree, mrca)
clr <- rep("black", nrow(tree$edge))
clr[which.edge(tree, desc)] <- "red"
plot(tree, edge.color = clr, show.tip.label = F, type = "fan")


Spatial data in R

For a spatial phylogenetic analysis, we need to define geographic areas of interest, within which organisms are considered part of the same community. These areas can be grid cells or can be irregular polygons, but we’ll focus on the former here as the latter are conceptually similar but more complex to work with.

We’ll be working with raster data and point data. The raster data comprise a spatial grid of presences and absences. Point data represent occurrence localities, which we’ll convert to grids of presences and absences before analysis.

Here’s a quick demo of loading and visualizing some spatial data.

# raster data, read in with the terra package
r <- rast("../../data/jetz_2012/bbs_occ.tif")
plot(r[[1:4]], col = "green")

plot(sum(r, na.rm = T))

# point data, read in as a data frame
p <- read_csv("../../data/mishler_2014/mishler_2014_acacia_points.csv")
## Rows: 132248 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Species, taxa
## dbl (4): Longitude, Latitude, x_metres_EPSG_3577_Albers_Equal_Area, y_metres...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# glimpse(p)
ggplot(p, aes(Longitude, Latitude)) + 
      geom_point()

# we can "rasterize" points by rounding the coordinates
# (there are fancier approaches but we'll keep it simple here)
pr <- p %>% mutate(x = round(Longitude), y = round(Latitude))
ggplot(pr, aes(x, y)) + 
      geom_raster()

# taking this a step farther, lets' convert our rounded coordinates 
# for one species into to a proper raster object
rp <- pr %>%
      filter(taxa == taxa[1000]) %>%
      dplyr::select(x, y) %>%
      distinct() %>%
      mutate(present = T) %>%
      rast()
plot(rp, col = "black")

Exercises:

  • Load the terrestrial mammal spatial data using rast()
  • Plot a map of the range of the gray wolf (Canis lupus).
  • Plot a map of canid species richness (all species in the genus Canis).
  • Plot a map of overall mammal species richness.
r <- rast("../../data/upham_2019/iucn_mammals.tif")
r %>% subset("Canis lupus") %>% plot()

r %>% subset(grep("Canis ", names(r))) %>% sum() %>% plot()

r %>% sum() %>% plot()


SPATIAL PHYLOGENETIC DATA

At a minimum, every spatial phylogenetic dataset consists of a phylogeny and a spatial dataset representing occurrence locations of the terminal taxa. We’ll use four different datasets as examples.

To keep things tidy, let’s start by loading the data for each one, and getting the pieces into a standardized format. We’ll package each dataset in a list with three pieces:

To avoid issues with some of the algorithms used later on, we need to do a few things to clean up the datasets, including removing taxa that aren’t present in both the spatial and phylogenetic datasets, and removing unoccupied sites from the community matrix. To do this we’ll use a function called clean_dataset() defined in functions.R.

Cailfornia vascular plants (from Thornhill et al. 2017):

# phylogeny
tree <- read.nexus("../../data/thornhill_2017/thornhill_2017_chronogram.nex")

# spatial data
occ <- read_csv("../../data/thornhill_2017/thornhill_2017_occ.csv")
## Rows: 1394170 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (1): otu
## dbl (2): x, y
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
glimpse(occ)
## Rows: 1,394,170
## Columns: 3
## $ x   <dbl> -121.8009, -121.0989, -122.2753, -121.6883, -121.2134, -119.5225, …
## $ y   <dbl> 37.39750, 40.07583, 40.74150, 37.21714, 40.02400, 37.44472, 33.533…
## $ otu <chr> "Delphinium", "Juncus", "Cyperus__Lipocarpha", "Cryptantha", "Cype…
# # compare species lists
all.equal(sort(tree$tip.label),
          sort(unique(occ$otu)))
## [1] TRUE
# summarize to lat-long grid
rnd <- function(x, y) round(x/y)*y # round x to the nearest y
occ <- occ %>%
      mutate(x = rnd(x, .5),
             y = rnd(y, .5)) %>%
      distinct()

# construct community matrix
comm <- occ %>%
      mutate(pres = 1) %>%
      spread(otu, pres, fill = 0)
xy <- select(comm, x, y)
comm <- select(comm, -x, -y) %>%
      as.matrix()

# package data
flora <- list(tree = tree,
              comm = comm,
              xy = xy) %>%
      clean_dataset()
## removing 221 sites with no occurrences

Australian acacia (from Mishler et al. 2014):

# phylogeny
tree <- read.tree("../../data/mishler_2014/mishler_2014_acacia.tre")

# spatial data
occ <- read_csv("../../data/mishler_2014/mishler_2014_acacia_points.csv")
## Rows: 132248 Columns: 6
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (2): Species, taxa
## dbl (4): Longitude, Latitude, x_metres_EPSG_3577_Albers_Equal_Area, y_metres...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# construct community matrix
occ <- occ %>%
      select(otu = taxa, x = Longitude, y = Latitude) %>%
      mutate(x = rnd(x, 1), # summarize to 0.5 degree lat-lon grid
             y = rnd(y, 1)) %>%
      distinct()
comm <- occ %>%
      mutate(pres = 1) %>%
      spread(otu, pres, fill = 0)
xy <- select(comm, x, y)
comm <- select(comm, -x, -y) %>%
      as.matrix()

# package data
acacia <- list(tree = tree,
               comm = comm,
               xy = xy) %>%
      clean_dataset()
## removing 797 sites with no occurrences

American birds (phylogeny from Jetz et al. 2012, spatial data from US Breeding Bird Survey):

# phylogeny (one of 10,000 posterior trees from Jetz 2012)
tree <- read.nexus("../../data/jetz_2012/jetz_2012_tree0001.nex")

# spatial data
r <- rast("../../data/jetz_2012/bbs_occ.tif")

# construct community matrix
comm <- values(r)
dim(comm)
## [1] 377 442
image(t(comm))

# spatial coordinates
xy <- crds(r[[1]], df = T, na.rm = F)

# package data
birds <- list(tree = tree,
              comm = comm,
              xy = xy) %>%
      clean_dataset()
## removing 242 sites with no occurrences

Mammals (phylogeny from Upham et al. 2019, spatial data from IUCN):

# phylogeny (one of 10,000 posterior trees from Upham 2019)
tree <- read.tree("../../data/upham_2019/upham_2019_tree0000.tre")

# spatial data
r <- rast("../../data/upham_2019/iucn_mammals.tif")

# construct community matrix
comm <- values(r)

# spatial coordinates
xy <- crds(r[[1]], df = T, na.rm = F)

# package data
mammals <- list(tree = tree,
                comm = comm,
                xy = xy) %>%
      clean_dataset()
## removing 734 tips from tree
## removing 446 taxa from community matrix
## removing 2598 sites with no occurrences

Exercises:

  • Plot the phylogeny for each of the datasets.
  • Plot the spatial coordinates for each of the datasets.
  • Identify the number of taxa and number of sites in each dataset.
mammals$tree %>% plot(show.tip.label = F)

mammals$xy %>% plot()

dim(mammals$comm)
## [1] 2598 5177

ALPHA DIVERSITY

Alpha phylogenetic diversity (and related metrics like phylogenetic endemism) represent the amount of evolutionary history present in a single community. We’ll look at functions from some existing packages that can calculate phylodiversity metrics, and reproduce these with some of our own calculations. We’ll also cover methods for significance testing, and explore the sensitivity of these tests to spatial scale. And we’ll take a look at how alternative evolutionary branch lengths influence results.

Pre-packaged functions

ds <- birds

# with phyloregion library
scomm <- ds$comm %>% Matrix(sparse = TRUE)
PD(scomm, ds$tree)
##   cell_1   cell_2   cell_3   cell_4   cell_5   cell_6   cell_7   cell_8 
## 3582.798 3778.919 4355.448 4312.056 3724.757 4320.454 3487.409 3529.021 
##   cell_9  cell_10  cell_11  cell_12  cell_13  cell_14  cell_15  cell_16 
## 3447.427 3419.588 3488.783 3588.469 4008.640 3696.477 4094.087 3811.433 
##  cell_17  cell_18  cell_19  cell_20  cell_21  cell_22  cell_23  cell_24 
## 3612.753 2590.026 2782.434 2385.522 2368.622 3627.100 3945.944 3931.083 
##  cell_25  cell_26  cell_27  cell_28  cell_29  cell_30  cell_31  cell_32 
## 3890.713 3331.068 3860.481 3590.577 3954.574 3948.325 3622.050 3517.718 
##  cell_33  cell_34  cell_35  cell_36  cell_37  cell_38  cell_39  cell_40 
## 3455.503 3687.643 3627.459 3894.832 4059.306 3866.207 3998.063 3652.541 
##  cell_41  cell_42  cell_43  cell_44  cell_45  cell_46  cell_47  cell_48 
## 3814.669 3610.482 3498.355 3519.732 3763.179 3905.573 4028.917 4138.715 
##  cell_49  cell_50  cell_51  cell_52  cell_53  cell_54  cell_55  cell_56 
## 4137.100 3728.158 3994.154 4252.226 4097.672 3939.590 4359.949 3402.430 
##  cell_57  cell_58  cell_59  cell_60  cell_61  cell_62  cell_63  cell_64 
## 3683.629 3767.627 3418.510 3660.785 3651.502 4020.410 3869.538 3785.533 
##  cell_65  cell_66  cell_67  cell_68  cell_69  cell_70  cell_71  cell_72 
## 3669.169 2827.882 2643.705 2877.501 3533.471 3870.439 3706.676 3722.440 
##  cell_73  cell_74  cell_75  cell_76  cell_77  cell_78  cell_79  cell_80 
## 3426.887 3804.753 4572.355 4429.111 3751.190 3880.011 4039.673 4363.419 
##  cell_81  cell_82  cell_83  cell_84  cell_85  cell_86  cell_87  cell_88 
## 3951.932 3878.154 4165.568 3871.093 3976.203 3591.010 3430.229 3321.550 
##  cell_89  cell_90  cell_91  cell_92  cell_93  cell_94  cell_95  cell_96 
## 3393.742 3060.813 3430.737 3654.015 3724.328 3467.054 3370.244 3401.804 
##  cell_97  cell_98  cell_99 cell_100 cell_101 cell_102 cell_103 cell_104 
## 3716.778 3698.316 3710.346 3718.237 3070.234 4019.254 4634.963 4472.986 
## cell_105 cell_106 cell_107 cell_108 cell_109 cell_110 cell_111 cell_112 
## 3642.810 3373.158 3784.174 4227.207 4147.057 4220.734 4151.165 3878.641 
## cell_113 cell_114 cell_115 cell_116 cell_117 cell_118 cell_119 cell_120 
## 3550.393 3085.749 3418.797 3564.853 3194.547 3272.234 3428.338 3436.437 
## cell_121 cell_122 cell_123 cell_124 cell_125 cell_126 cell_127 cell_128 
## 3061.708 3131.114 3212.330 3462.690 3467.245 3814.418 3584.872 2618.477 
## cell_129 cell_130 cell_131 cell_132 cell_133 cell_134 cell_135 cell_136 
## 3392.011 4436.385 4745.206 3523.526 3417.811 3578.737 4255.684 3914.282 
## cell_137 cell_138 cell_139 cell_140 cell_141 cell_142 cell_143 cell_144 
## 4240.924 4344.170 4166.059 3559.037 3544.248 4062.733 3311.566 3478.943 
## cell_145 cell_146 cell_147 cell_148 cell_149 cell_150 cell_151 cell_152 
## 3205.155 3462.813 3526.216 3311.393 2976.729 2900.984 3297.404 3167.508 
## cell_153 cell_154 cell_155 cell_156 cell_157 cell_158 cell_159 cell_160 
## 3906.271 3764.395 4224.789 4462.119 3281.642 3052.589 3343.311 2659.387 
## cell_161 cell_162 cell_163 cell_164 cell_165 cell_166 cell_167 cell_168 
## 3311.178 4237.474 3848.531 3315.250 3692.636 3543.209 3303.631 3288.513 
## cell_169 cell_170 cell_171 cell_172 cell_173 cell_174 cell_175 cell_176 
## 3177.494 3376.268 3261.661 3252.612 3275.281 3363.141 3140.545 3421.232 
## cell_177 cell_178 cell_179 cell_180 cell_181 cell_182 cell_183 cell_184 
## 3561.333 3325.336 4315.780 4102.393 3150.967 4092.881 4302.440 3293.989 
## cell_185 cell_186 cell_187 cell_188 cell_189 cell_190 cell_191 cell_192 
## 4099.675 3108.686 3229.429 3500.531 3682.134 3409.398 3657.392 3579.778 
## cell_193 cell_194 cell_195 cell_196 cell_197 cell_198 cell_199 cell_200 
## 3605.348 3484.089 3529.832 3409.855 3404.296 3531.209 3737.846 3355.301 
## cell_201 cell_202 cell_203 cell_204 cell_205 cell_206 cell_207 cell_208 
## 2280.667 3641.156 3023.014 3442.615 4088.806 3506.021 3845.148 2944.231 
## cell_209 cell_210 cell_211 cell_212 cell_213 cell_214 cell_215 cell_216 
## 2889.317 3297.106 3285.397 3741.039 3540.055 3673.718 3468.004 3366.721 
## cell_217 cell_218 cell_219 cell_220 cell_221 cell_222 cell_223 cell_224 
## 3555.077 3503.021 3688.419 3532.081 3582.126 2974.381 3442.856 3882.064 
## cell_225 cell_226 cell_227 cell_228 cell_229 cell_230 cell_231 cell_232 
## 4146.162 4184.051 3906.099 4019.940 3700.945 3966.544 3752.479 3519.111 
## cell_233 cell_234 cell_235 cell_236 cell_237 cell_238 cell_239 cell_240 
## 3489.316 4382.977 4118.758 4066.206 3878.009 2474.686 4377.796 3679.241 
## cell_241 cell_242 
## 4177.109 1966.891
pr <- ds$xy %>%
      mutate(pd = PD(scomm, ds$tree),
             pe = phylo_endemism(scomm, ds$tree),
             we = weighted_endemism(scomm))
pr %>% ggplot(aes(x, y, fill = pd)) + geom_raster()

pr %>% carto("pd")

pr %>% carto("pe")

# with canaper library
# ?cpr_rand_test
cpr <- cpr_rand_test(ds$comm, ds$tree, 
                     null_model = "curveball",
                     n_reps = 1) %>%
      bind_cols(ds$xy)
# glimpse(cpr)
# names(cpr)
cpr %>% carto("pd_obs")

cpr %>% carto("pe_obs")

# plot a local community on the phylogeny
par(mfrow = c(1, 2))
for(fun in list(min, max)){
      cell <- cpr %>% filter(pd_obs == fun(pd_obs)) %>% rownames()
      spp <- colnames(ds$comm)[ds$comm[cell,] == 1]
      phyplot(ds$tree, connect = spp, 
              type = "fan", show.tip.label = F,
              main = paste("taxa found in", cell))
}

Manual approach

Canaper and phyloregion are pretty great, but how are these spatial biodiversity metrics actually being calculated under the hood? Let’s reproduce some of these results with a slightly more manual approach.

# build community matrix including not just tips as above,
# but also every clade at every level
pcomm <- build_clade_ranges(ds$tree, ds$comm) # from functions.r
# dim(ds$comm)
# dim(pcomm)
image(t(pcomm))

# clade-level metrics
range_size <- colSums(pcomm)
branch_length <- ds$tree$edge.length / sum(ds$tree$edge.length)
phyplot(ds$tree, value = branch_length,
        show.tip.label = F, type = "fan", main = "branch length")

phyplot(ds$tree, value = log(1 / range_size),
        show.tip.label = F, type = "fan", main = "log endemism")

# PD and PE are just weighted summaries across this pcomm matrix:
pd <- pcomm %>%
      apply(1, function(x) x * branch_length) %>%
      colSums()
ds$xy %>% mutate(pd = pd) %>% carto("pd")

pe <- pcomm %>%
      apply(1, function(x) x * branch_length / range_size) %>%
      colSums()
ds$xy %>% mutate(pe = pe) %>% carto("pe")

Exercises:

  • Compute these metrics with one of the other datasets.
  • Species richness can be considered PD on a star phylogeny. Prove it.
  • Manual calculations make it clear that these diversity metrics are simply weighted summaries of the pcomm matrix, which invites any number of additional weighting schemes. Make a version weighted by: taxon endangerment; site endangerment; or taxon-site habitat importance (using fake simulated data for these variables).
# SR == PD on star phylogeny
ds <- birds
star <- ds$tree
tip <- ! star$edge[,2] %in% star$edge[,1]
star$edge.length <- rep(0, length(star$edge.length))
star$edge.length[tip] <- 1
plot(star, type = "fan", show.tip.label = F)

sr <- rowSums(ds$comm)
pd <- PD(scomm, star)
plot(sr, pd)

# taxon & site endangerment
# (this would represent locations of treatened sites with high concentrations of range-restricted, endangered taxa)
endangerment <- runif(ncol(pcomm)) # fake random data
site_threat <- runif(nrow(pcomm)) # fake random data
tpe <- pcomm %>%
      apply(1, function(x) x * branch_length / range_size * endangerment) %>%
      apply(1, function(x) x * site_threat) %>%
      rowSums()
ds$xy %>% mutate(tpe = tpe) %>% carto("tpe")

Significance testing

To get a sense for how noteworthy an observed high or low diversity metric is, it can be helpful to compare the observed value to a null distribution. Let’s demonstrate this using the acacia dataset. The cpr_rand_test() function randomly permutes the community matrix many times, calculating the diversity of each location each time to derive a distribution of null diversity values for each location.

CANAPE (Mishler et al. 2014) is a special application of this randomization approach, used to identify areas of neo-endemism and paleo-endemism.

ds <- acacia

# note that this is a phylogram, not a chronogram
# ds$tree %>% plot(show.tip.label = F)

# phylodiversity randomizations
cpr <- cpr_rand_test(ds$comm, ds$tree, 
                     null_model = "curveball",
                     n_reps = 101, n_iterations = 10000) %>%
      bind_cols(ds$xy)

# raw PD patterns here are very different from randomized quantiles
cpr %>% carto("pd_obs")

cpr %>% carto("pd_obs_p_upper")

cpr %>%
      ggplot(aes(pd_obs, pd_obs_p_upper)) +
      geom_point()

cpr %>% carto("rpd_obs_p_upper")

## CANAPE ##

# classify canape categories
end <- cpr %>% cpr_classify_endem() 

# map
end %>%
      carto("endem_type") +
      scale_fill_manual(
            values = c("mediumpurple1", "red", "gray80", "blue", "darkorchid4"),
            breaks = c("mixed", "neo", "not significant", "paleo", "super"))

# scatterplot of canape components
end %>%
      ggplot(aes(pe_alt_obs, pe_obs,
                 color = endem_type)) +
      geom_point() +
      scale_color_manual(
            values = c("mediumpurple1", "red", "gray80", "blue", "darkorchid4"),
            breaks = c("mixed", "neo", "not significant", "paleo", "super"))

Null models can be constructed in many different ways, and these differences can have a major effect on results. We won’t get too far into the weeds here, but let’s illustrate how a few different null models affect PD significance values for this dataset:

# ?vegan::commsim
curveball <- cpr %>% 
      mutate(model = "curveball")
r0 <- cpr_rand_test(ds$comm, ds$tree, null_model = "r0", metrics = "pd", 
                    n_reps = 101) %>%
      bind_cols(ds$xy) %>%
      mutate(model = "r0")
c0 <- cpr_rand_test(ds$comm, ds$tree, null_model = "c0", metrics = "pd", 
                    n_reps = 101) %>%
      bind_cols(ds$xy) %>%
      mutate(model = "c0")

bind_rows(curveball, r0, c0) %>%
      mutate(pd_sig = case_when(pd_obs_p_upper > .95 ~ "high",
                                pd_obs_p_upper < .05 ~ "low",
                                TRUE ~ "no")) %>%
      carto("pd_sig") + 
      facet_wrap(~model, nrow = 1) +
      scale_fill_manual(values = c("red", "blue", "gray"))

Before moving on from randomizations, it’s important to note that the size of the geographic region in an analysis can strongly influence the results, because it alters the species pool to which a site is compared during randomization. Let’s look at three different extents:

# a function to subset the sites in a dataset
filter_sites <- function(ds, selection){
      dsf <- ds
      dsf$comm <- dsf$comm[selection, ]
      dsf$xy <- dsf$xy[selection, ]
      dsf<- clean_dataset(dsf)
      return(dsf)
}

# filter mammals data to three nested regions of different sizes
americas <- mammals %>% 
      filter_sites(mammals$xy$x < -34)
## removing 3415 taxa with no occurrences
## removing 3415 tips from tree
## removing 914 sites with no occurrences
s_america <- mammals %>% 
      filter_sites(between(mammals$xy$x, -82, -34) 
                   & mammals$xy$y < 13)
## removing 4027 taxa with no occurrences
## removing 4027 tips from tree
## removing 227 sites with no occurrences
amazon <- mammals %>% 
      filter_sites(between(mammals$xy$x, -75, -50) 
                   & between(mammals$xy$y, -15, 5))
## removing 4562 taxa with no occurrences
## removing 4562 tips from tree
## removing 56 sites with no occurrences
# compute pd significance for each dataset
cpr_americas <- cpr_rand_test(
      americas$comm, americas$tree, 
      null_model = "curveball", metrics = "pd", n_reps = 21) %>%
      bind_cols(americas$xy) %>%
      mutate(region = "americas")
cpr_s_america <- cpr_rand_test(
      s_america$comm, s_america$tree, 
      null_model = "curveball", metrics = "pd", n_reps = 21) %>%
      bind_cols(s_america$xy) %>%
      mutate(region = "south america")
cpr_amazon <- cpr_rand_test(
      amazon$comm, amazon$tree, 
      null_model = "curveball", metrics = "pd", n_reps = 21) %>%
      bind_cols(amazon$xy) %>%
      mutate(region = "amazon")
cpr <- bind_rows(cpr_americas, cpr_s_america, cpr_amazon) %>%
      mutate(region = factor(region, levels = unique(region))) %>%
      mutate(pd_sig = case_when(pd_obs_p_upper > .95 ~ "high",
                                pd_obs_p_upper < .05 ~ "low",
                                TRUE ~ "no"))

# maps
cpr %>%
      carto("pd_sig") +
      facet_wrap(~ region, nrow = 1) +
      scale_fill_manual(values = c("red", "blue", "gray"))

cpr %>%
      filter(x %in% x[region == "amazon"],
             y %in% y[region == "amazon"]) %>%
      ggplot(aes(x, y, fill = pd_sig)) +
      geom_raster() +
      facet_wrap(~ region, nrow = 1) +
      theme_void() + theme(legend.position = "top") + coord_fixed() +
      scale_fill_manual(values = c("red", "blue", "gray"))

Facets of phylogenetic diversity

Phylogenetic branch lengths can represent different macroevolutionary variables, each with different implications for biogeography and conservation. We’ve indirectly met some different “phylogenetic diversity facets” already, through RPD (which is computed in canaper). Let’s quickly explore them in a bit more detail.

ds <- flora

# function to normalize branch lengths
scale_edges <- function(tree){
      tree$edge.length <- tree$edge.length / sum(tree$edge.length)
      return(tree)
}

# load trees
chronogram <- read.nexus("../../data/thornhill_2017/thornhill_2017_chronogram.nex") %>%
      scale_edges()
phylogram <- read.nexus("../../data/thornhill_2017/thornhill_2017_phylogram.nex") %>%
      scale_edges()

# make cladogram
cladogram <- phylogram
cladogram$edge.length <- rep(1, length(cladogram$edge.length))
cladogram <- scale_edges(cladogram)

# plot trees
par(mfrow = c(1, 3))
phyplot(chronogram, chronogram$edge.length,
        type = "fan", show.tip.label = F, main = "time")
phyplot(phylogram, phylogram$edge.length,
        type = "fan", show.tip.label = F, main = "divergence")
phyplot(cladogram, cladogram$edge.length,
        type = "fan", show.tip.label = F, main = "cladogenesis")

dev.off()
## null device 
##           1
# format as sparse matrix for phyloregion package
scomm <- ds$comm %>% Matrix(sparse = TRUE)

# PD for the three trees
div <- ds$xy %>%
      mutate(divergence = PD(scomm, phylogram),
             time = PD(scomm, chronogram),
             diversification = PD(scomm, cladogram))

# visualize
div %>%
      pivot_longer(c(divergence, time, diversification), 
                   "facet", "value") %>%
      group_by(facet) %>%
      mutate(rank = rank(value)) %>%
      carto("rank") + 
      facet_wrap(~facet, nrow = 1)
div %>%
      mutate(rpd = time / diversification) %>% 
      carto("rpd")
div %>%
      mutate(dynamism = divergence / time) %>%
      carto("dynamism")

Exercises:

  • Verify that rpd as returned by the cpr_rand_test() function using a chronogram is identical to the ratio of chronogram PD to cladogram PD calculated immediately above.
  • Make a map of the ratio of phylogram PD to cladogram PD, which measures the relative contribution of anagenesis versus cladogenesis to a community’s history. Where is this highest and lowest?
# RPD from cpr_rand_test() == chronogram PD / cladogram PD
cpr <- cpr_rand_test(ds$comm, ds$tree, null_model = "curveball", n_reps = 1)
pd_chron <- PD(scomm, chronogram)
pd_clad <- PD(scomm, cladogram)
plot(div$time/div$diversification, cpr$rpd_obs)

# phylogram PD / cladogram PD
div %>%
      mutate(genesis = divergence / diversification) %>%
      carto("genesis")

# phylogram PE / cladogram PE
ds$xy %>%
      mutate(divergence = phylo_endemism(scomm, phylogram),
             time = phylo_endemism(scomm, chronogram),
             diversification = phylo_endemism(scomm, cladogram)) %>%
      mutate(genesis = divergence / diversification) %>%
      carto("genesis")


BETA DIVERSITY

Analyses of beta diversity focus on compositional differences among local communities. Phylogenetic versions of these analyses represent differences between communities as the difference in the portion of the phylogenetic tree that’s found in each site.

Turnover

We’ll begin by creating a pairwise turnover matrix with a row and column for each site. We can compare this to species-based turnover and geographic distance. Further exploration of these turnover patterns could use a modeling approach like GDM (Ferrier 2007) to understand how factors like environmental differences explain phylogenetic turnover, but we will not cover that here.

ds <- mammals

# construct dissimlarity matrices
phy_beta <- ds$comm %>% # phylo sorensen's
      Matrix(sparse = TRUE) %>%
      phylobeta(ds$tree) %>%
      pluck("phylo.beta.sor") 
str(phy_beta)
##  'dist' Named num [1:3373503] 0.123 0.23 0.23 0.23 0.23 ...
##  - attr(*, "Labels")= chr [1:2598] "cell_1" "cell_2" "cell_3" "cell_4" ...
##  - attr(*, "Size")= int 2598
##  - attr(*, "Diag")= logi FALSE
##  - attr(*, "Upper")= logi FALSE
sp_beta <- ds$comm %>% # species sorensen's
      betapart::beta.pair() %>%
      pluck("beta.sor")
geo_dist <- ds$xy %>% # geographic distance
      geosphere::distm()

# plot geographic distance against phylogenetic turnover
ss <- sample(length(geo_dist), 10000)
tibble(distance = geo_dist[ss],
       phylo_sor = as.matrix(phy_beta)[ss],
       species_sor = as.matrix(sp_beta)[ss]) %>%
      pivot_longer(-distance, "stat", "value") %>%
      ggplot(aes(distance, value, color = stat)) +
      geom_point(size = .5) +
      geom_smooth(aes(group = stat), color = "black")
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Regionalization

We can also use turnover data to make maps of phylogenetic similarity. Let’s look at two approaches: a discrete classification of phylogenetic regions using a cluster analysis from the phyloregions package, and a continuous ordination map that plots similar sites in similar colors:

# phylogenetic regionalization (spatial cluster analysis)
k <- 25
phyloregion(phy_beta, k = k) %>%
      pluck("membership") %>%
      bind_cols(ds$xy) %>%
      mutate(cluster = factor(cluster)) %>%
      carto("cluster") +
      scale_fill_manual(values = unname(pals::alphabet(k)))

# nmds of regions
    reg <- phyloregion(phy_beta, k = k)
    metaMDS(reg$region.dist, k = 2)$points %>%
      as.data.frame() %>%
      ggplot(aes(MDS1, MDS2, color = factor(1:k), label = 1:k)) +
      geom_label(fontface = "bold") +
      scale_color_manual(values = unname(pals::alphabet(k))) +
      theme(legend.position = "none")
## Run 0 stress 0.2018687 
## Run 1 stress 0.1874266 
## ... New best solution
## ... Procrustes: rmse 0.1235785  max resid 0.3792318 
## Run 2 stress 0.2208355 
## Run 3 stress 0.1975499 
## Run 4 stress 0.1874266 
## ... New best solution
## ... Procrustes: rmse 6.565778e-05  max resid 0.0001897589 
## ... Similar to previous best
## Run 5 stress 0.2001864 
## Run 6 stress 0.2222201 
## Run 7 stress 0.2243932 
## Run 8 stress 0.2464383 
## Run 9 stress 0.1997677 
## Run 10 stress 0.2018926 
## Run 11 stress 0.2296424 
## Run 12 stress 0.1997677 
## Run 13 stress 0.1975499 
## Run 14 stress 0.1893562 
## Run 15 stress 0.1991909 
## Run 16 stress 0.1893562 
## Run 17 stress 0.2633435 
## Run 18 stress 0.1975499 
## Run 19 stress 0.1874266 
## ... Procrustes: rmse 1.75014e-06  max resid 5.302682e-06 
## ... Similar to previous best
## Run 20 stress 0.2166999 
## *** Solution reached

    # unpacking and visualizing regional hclust
    library(dendextend)
## Registered S3 method overwritten by 'dendextend':
##   method     from 
##   rev.hclust vegan
## 
## ---------------------
## Welcome to dendextend version 1.16.0
## Type citation('dendextend') for how to cite the package.
## 
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
## 
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags: 
##   https://stackoverflow.com/questions/tagged/dendextend
## 
##  To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
## ---------------------
## 
## Attaching package: 'dendextend'
## The following object is masked from 'package:permute':
## 
##     shuffle
## The following object is masked from 'package:raster':
## 
##     rotate
## The following object is masked from 'package:terra':
## 
##     rotate
## The following objects are masked from 'package:ape':
## 
##     ladderize, rotate
## The following object is masked from 'package:stats':
## 
##     cutree
    hc <- phy_beta %>% hclust(method = "average")
    clust <- phyloregion(phy_beta, k = k) %>%
          pluck("membership")
    clrs <- unname(pals::alphabet(k))[clust$cluster]
    dend <- as.dendrogram(hc)
    o <- order.dendrogram(dend)
    dend <- assign_values_to_leaves_edgePar(
          dend, value = clrs[o], edgePar = "col")
    plot(dend, leaflab = "none")

    circlize_dendrogram(dend, labels = F, 
                        dend_track_height = .9, 
                        labels_track_height = 0)
## Loading required namespace: circlize

# RGB visualization
# ord <- vegan::metaMDS(phy_beta, k = 3, trymax = 50)$points
# saveRDS(ord$points, "../../data/cache/mammals_ord.rds")
ord <- readRDS("../../data/cache/mammals_ord.rds")
ord %>%
      as.data.frame() %>%
      mutate_all(function(x) rank(x)/max(rank(x))) %>%
      bind_cols(ds$xy) %>%
      mutate(color = rgb(MDS1, MDS2, MDS3)) %>%
      carto("color") +
      scale_fill_identity()

Exercises:

  • The above demonstration used Sorensen’s distance as a measure of phylogenetic community difference. Try things with Jaccard’s distance. How do the results differ?
  • Experiment with different numbers of clusters and different cluster methods in the phyloregion() function, to explore how robust regional boundaries are to these parameter choices.

CONSERVATION PRIORITIZATION

Phylogenetic diversity is useful for basic research, but from the beginning it has also been proposed as an applied conservation tool. Let’s use it to run a conservation prioritization analysis. We’ll focus on the California vascular flora dataset for this exercise, which has been used for conservation prioritization in the past (Kling et al. 2018).

Optimal reserve design is a major computational challenge, due to the extremely large number of sets of sites that could comprise a protected area network. Our computational workhorse will be the prioritizr library, a very flexible and powerful conservation planning toolkit that uses linear solvers to identify optimal conservation solutions.

First we need to get our data formatted as raster objects. We’ll also load an additional data layer of current protected areas, derived from Kling et al. (2018).

ds <- flora

# conservation features
features <- ds$xy %>%
      bind_cols(ds$comm) %>%
      rasterFromXYZ()

# spatial planning units
units <- features[[1]] %>% 
      reclassify(c(-Inf, Inf, 1))

# current protected areas
protected <- raster("../../data/thornhill_2017/kling_2018_protection_status.tif")
# plot(protected)
protected <- protected > .66
prot_area <- sum(protected[], na.rm = T)

# harmonize spatial coverage
overlap <- protected + units
protected <- protected %>% crop(overlap) %>% mask(overlap)
units <- units %>% crop(overlap) %>% mask(overlap)
features <- features %>% crop(overlap) %>% mask(overlap)
crs(protected) <- crs(units) <- crs(features)

The prioritizr library has some native functionality for incorporating phylogenies into conservation prioritization, so let’s start there.

# compare existing protected area network to 
# an optimal network of the same size
prob <- problem(units, features) %>%
      add_max_phylo_div_objective(budget = prot_area,
                                  tree = ds$tree) %>%
      add_relative_targets(0.1) %>%
      add_binary_decisions() %>%
      add_rsymphony_solver(gap = 0)
soln <- solve(prob)
stack(soln, protected) %>% 
      setNames(c("optimal", "existing")) %>%
      plot(col = c("gray", "orange"))

# identify the highest-priority sites to add to existing protected areas
prob <- prob %>%
      add_max_phylo_div_objective(budget = prot_area + 10,
                                  tree = ds$tree) %>%
      add_locked_in_constraints(locked_in = protected)
## Warning in res(x, ...): overwriting previously defined objective
soln <- solve(prob)
plot(soln + protected, 
     col = c("gray", "orange", "darkgreen"))

While it incorporates phylogeny, this method is still “tip-centric” in the sense that its objective is to maximize the phylogenetic diversity of the tips protected across the reserve network – the targets are still the tips, rather than taxa at all levels. A subtly different alternative that’s arguably more consistent with spatial phylogenetic metrics like PE is to build out a phylogenetic community matrix, and weight each taxon by its branch length during prioritization. We can see this yields a somewhat different result:

# build feature set representing all clades, not just tips
pcomm <- build_clade_ranges(ds$tree, ds$comm)
pfeatures <- ds$xy %>%
      bind_cols(pcomm) %>%
      rasterFromXYZ()
## New names:
## • `` -> `...3`
## • `` -> `...4`
## • `` -> `...5`
## • `` -> `...6`
## • `` -> `...7`
## • `` -> `...8`
## • `` -> `...9`
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
## • `` -> `...13`
## • `` -> `...14`
## • `` -> `...15`
## • `` -> `...16`
## • `` -> `...17`
## • `` -> `...18`
## • `` -> `...19`
## • `` -> `...20`
## • `` -> `...21`
## • `` -> `...22`
## • `` -> `...23`
## • `` -> `...24`
## • `` -> `...25`
## • `` -> `...26`
## • `` -> `...27`
## • `` -> `...28`
## • `` -> `...29`
## • `` -> `...30`
## • `` -> `...31`
## • `` -> `...32`
## • `` -> `...33`
## • `` -> `...34`
## • `` -> `...35`
## • `` -> `...36`
## • `` -> `...37`
## • `` -> `...38`
## • `` -> `...39`
## • `` -> `...40`
## • `` -> `...41`
## • `` -> `...42`
## • `` -> `...43`
## • `` -> `...44`
## • `` -> `...45`
## • `` -> `...46`
## • `` -> `...47`
## • `` -> `...48`
## • `` -> `...49`
## • `` -> `...50`
## • `` -> `...51`
## • `` -> `...52`
## • `` -> `...53`
## • `` -> `...54`
## • `` -> `...55`
## • `` -> `...56`
## • `` -> `...57`
## • `` -> `...58`
## • `` -> `...59`
## • `` -> `...60`
## • `` -> `...61`
## • `` -> `...62`
## • `` -> `...63`
## • `` -> `...64`
## • `` -> `...65`
## • `` -> `...66`
## • `` -> `...67`
## • `` -> `...68`
## • `` -> `...69`
## • `` -> `...70`
## • `` -> `...71`
## • `` -> `...72`
## • `` -> `...73`
## • `` -> `...74`
## • `` -> `...75`
## • `` -> `...76`
## • `` -> `...77`
## • `` -> `...78`
## • `` -> `...79`
## • `` -> `...80`
## • `` -> `...81`
## • `` -> `...82`
## • `` -> `...83`
## • `` -> `...84`
## • `` -> `...85`
## • `` -> `...86`
## • `` -> `...87`
## • `` -> `...88`
## • `` -> `...89`
## • `` -> `...90`
## • `` -> `...91`
## • `` -> `...92`
## • `` -> `...93`
## • `` -> `...94`
## • `` -> `...95`
## • `` -> `...96`
## • `` -> `...97`
## • `` -> `...98`
## • `` -> `...99`
## • `` -> `...100`
## • `` -> `...101`
## • `` -> `...102`
## • `` -> `...103`
## • `` -> `...104`
## • `` -> `...105`
## • `` -> `...106`
## • `` -> `...107`
## • `` -> `...108`
## • `` -> `...109`
## • `` -> `...110`
## • `` -> `...111`
## • `` -> `...112`
## • `` -> `...113`
## • `` -> `...114`
## • `` -> `...115`
## • `` -> `...116`
## • `` -> `...117`
## • `` -> `...118`
## • `` -> `...119`
## • `` -> `...120`
## • `` -> `...121`
## • `` -> `...122`
## • `` -> `...123`
## • `` -> `...124`
## • `` -> `...125`
## • `` -> `...126`
## • `` -> `...127`
## • `` -> `...128`
## • `` -> `...129`
## • `` -> `...130`
## • `` -> `...131`
## • `` -> `...132`
## • `` -> `...133`
## • `` -> `...134`
## • `` -> `...135`
## • `` -> `...136`
## • `` -> `...137`
## • `` -> `...138`
## • `` -> `...139`
## • `` -> `...140`
## • `` -> `...141`
## • `` -> `...142`
## • `` -> `...143`
## • `` -> `...144`
## • `` -> `...145`
## • `` -> `...146`
## • `` -> `...147`
## • `` -> `...148`
## • `` -> `...149`
## • `` -> `...150`
## • `` -> `...151`
## • `` -> `...152`
## • `` -> `...153`
## • `` -> `...154`
## • `` -> `...155`
## • `` -> `...156`
## • `` -> `...157`
## • `` -> `...158`
## • `` -> `...159`
## • `` -> `...160`
## • `` -> `...161`
## • `` -> `...162`
## • `` -> `...163`
## • `` -> `...164`
## • `` -> `...165`
## • `` -> `...166`
## • `` -> `...167`
## • `` -> `...168`
## • `` -> `...169`
## • `` -> `...170`
## • `` -> `...171`
## • `` -> `...172`
## • `` -> `...173`
## • `` -> `...174`
## • `` -> `...175`
## • `` -> `...176`
## • `` -> `...177`
## • `` -> `...178`
## • `` -> `...179`
## • `` -> `...180`
## • `` -> `...181`
## • `` -> `...182`
## • `` -> `...183`
## • `` -> `...184`
## • `` -> `...185`
## • `` -> `...186`
## • `` -> `...187`
## • `` -> `...188`
## • `` -> `...189`
## • `` -> `...190`
## • `` -> `...191`
## • `` -> `...192`
## • `` -> `...193`
## • `` -> `...194`
## • `` -> `...195`
## • `` -> `...196`
## • `` -> `...197`
## • `` -> `...198`
## • `` -> `...199`
## • `` -> `...200`
## • `` -> `...201`
## • `` -> `...202`
## • `` -> `...203`
## • `` -> `...204`
## • `` -> `...205`
## • `` -> `...206`
## • `` -> `...207`
## • `` -> `...208`
## • `` -> `...209`
## • `` -> `...210`
## • `` -> `...211`
## • `` -> `...212`
## • `` -> `...213`
## • `` -> `...214`
## • `` -> `...215`
## • `` -> `...216`
## • `` -> `...217`
## • `` -> `...218`
## • `` -> `...219`
## • `` -> `...220`
## • `` -> `...221`
## • `` -> `...222`
## • `` -> `...223`
## • `` -> `...224`
## • `` -> `...225`
## • `` -> `...226`
## • `` -> `...227`
## • `` -> `...228`
## • `` -> `...229`
## • `` -> `...230`
## • `` -> `...231`
## • `` -> `...232`
## • `` -> `...233`
## • `` -> `...234`
## • `` -> `...235`
## • `` -> `...236`
## • `` -> `...237`
## • `` -> `...238`
## • `` -> `...239`
## • `` -> `...240`
## • `` -> `...241`
## • `` -> `...242`
## • `` -> `...243`
## • `` -> `...244`
## • `` -> `...245`
## • `` -> `...246`
## • `` -> `...247`
## • `` -> `...248`
## • `` -> `...249`
## • `` -> `...250`
## • `` -> `...251`
## • `` -> `...252`
## • `` -> `...253`
## • `` -> `...254`
## • `` -> `...255`
## • `` -> `...256`
## • `` -> `...257`
## • `` -> `...258`
## • `` -> `...259`
## • `` -> `...260`
## • `` -> `...261`
## • `` -> `...262`
## • `` -> `...263`
## • `` -> `...264`
## • `` -> `...265`
## • `` -> `...266`
## • `` -> `...267`
## • `` -> `...268`
## • `` -> `...269`
## • `` -> `...270`
## • `` -> `...271`
## • `` -> `...272`
## • `` -> `...273`
## • `` -> `...274`
## • `` -> `...275`
## • `` -> `...276`
## • `` -> `...277`
## • `` -> `...278`
## • `` -> `...279`
## • `` -> `...280`
## • `` -> `...281`
## • `` -> `...282`
## • `` -> `...283`
## • `` -> `...284`
## • `` -> `...285`
## • `` -> `...286`
## • `` -> `...287`
## • `` -> `...288`
## • `` -> `...289`
## • `` -> `...290`
## • `` -> `...291`
## • `` -> `...292`
## • `` -> `...293`
## • `` -> `...294`
## • `` -> `...295`
## • `` -> `...296`
## • `` -> `...297`
## • `` -> `...298`
## • `` -> `...299`
## • `` -> `...300`
## • `` -> `...301`
## • `` -> `...302`
## • `` -> `...303`
## • `` -> `...304`
## • `` -> `...305`
## • `` -> `...306`
## • `` -> `...307`
## • `` -> `...308`
## • `` -> `...309`
## • `` -> `...310`
## • `` -> `...311`
## • `` -> `...312`
## • `` -> `...313`
## • `` -> `...314`
## • `` -> `...315`
## • `` -> `...316`
## • `` -> `...317`
## • `` -> `...318`
## • `` -> `...319`
## • `` -> `...320`
## • `` -> `...321`
## • `` -> `...322`
## • `` -> `...323`
## • `` -> `...324`
## • `` -> `...325`
## • `` -> `...326`
## • `` -> `...327`
## • `` -> `...328`
## • `` -> `...329`
## • `` -> `...330`
## • `` -> `...331`
## • `` -> `...332`
## • `` -> `...333`
## • `` -> `...334`
## • `` -> `...335`
## • `` -> `...336`
## • `` -> `...337`
## • `` -> `...338`
## • `` -> `...339`
## • `` -> `...340`
## • `` -> `...341`
## • `` -> `...342`
## • `` -> `...343`
## • `` -> `...344`
## • `` -> `...345`
## • `` -> `...346`
## • `` -> `...347`
## • `` -> `...348`
## • `` -> `...349`
## • `` -> `...350`
## • `` -> `...351`
## • `` -> `...352`
## • `` -> `...353`
## • `` -> `...354`
## • `` -> `...355`
## • `` -> `...356`
## • `` -> `...357`
## • `` -> `...358`
## • `` -> `...359`
## • `` -> `...360`
## • `` -> `...361`
## • `` -> `...362`
## • `` -> `...363`
## • `` -> `...364`
## • `` -> `...365`
## • `` -> `...366`
## • `` -> `...367`
## • `` -> `...368`
## • `` -> `...369`
## • `` -> `...370`
## • `` -> `...371`
## • `` -> `...372`
## • `` -> `...373`
## • `` -> `...374`
## • `` -> `...375`
## • `` -> `...376`
## • `` -> `...377`
## • `` -> `...378`
## • `` -> `...379`
## • `` -> `...380`
## • `` -> `...381`
## • `` -> `...382`
## • `` -> `...383`
## • `` -> `...384`
## • `` -> `...385`
## • `` -> `...386`
## • `` -> `...387`
## • `` -> `...388`
## • `` -> `...389`
## • `` -> `...390`
## • `` -> `...391`
## • `` -> `...392`
## • `` -> `...393`
## • `` -> `...394`
## • `` -> `...395`
## • `` -> `...396`
## • `` -> `...397`
## • `` -> `...398`
## • `` -> `...399`
## • `` -> `...400`
## • `` -> `...401`
## • `` -> `...402`
## • `` -> `...403`
## • `` -> `...404`
## • `` -> `...405`
## • `` -> `...406`
## • `` -> `...407`
## • `` -> `...408`
## • `` -> `...409`
## • `` -> `...410`
## • `` -> `...411`
## • `` -> `...412`
## • `` -> `...413`
## • `` -> `...414`
## • `` -> `...415`
## • `` -> `...416`
## • `` -> `...417`
## • `` -> `...418`
## • `` -> `...419`
## • `` -> `...420`
## • `` -> `...421`
## • `` -> `...422`
## • `` -> `...423`
## • `` -> `...424`
## • `` -> `...425`
## • `` -> `...426`
## • `` -> `...427`
## • `` -> `...428`
## • `` -> `...429`
## • `` -> `...430`
## • `` -> `...431`
## • `` -> `...432`
## • `` -> `...433`
## • `` -> `...434`
## • `` -> `...435`
## • `` -> `...436`
## • `` -> `...437`
## • `` -> `...438`
## • `` -> `...439`
## • `` -> `...440`
## • `` -> `...441`
## • `` -> `...442`
## • `` -> `...443`
## • `` -> `...444`
## • `` -> `...445`
## • `` -> `...446`
## • `` -> `...447`
## • `` -> `...448`
## • `` -> `...449`
## • `` -> `...450`
## • `` -> `...451`
## • `` -> `...452`
## • `` -> `...453`
## • `` -> `...454`
## • `` -> `...455`
## • `` -> `...456`
## • `` -> `...457`
## • `` -> `...458`
## • `` -> `...459`
## • `` -> `...460`
## • `` -> `...461`
## • `` -> `...462`
## • `` -> `...463`
## • `` -> `...464`
## • `` -> `...465`
## • `` -> `...466`
## • `` -> `...467`
## • `` -> `...468`
## • `` -> `...469`
## • `` -> `...470`
## • `` -> `...471`
## • `` -> `...472`
## • `` -> `...473`
## • `` -> `...474`
## • `` -> `...475`
## • `` -> `...476`
## • `` -> `...477`
## • `` -> `...478`
## • `` -> `...479`
## • `` -> `...480`
## • `` -> `...481`
## • `` -> `...482`
## • `` -> `...483`
## • `` -> `...484`
## • `` -> `...485`
## • `` -> `...486`
## • `` -> `...487`
## • `` -> `...488`
## • `` -> `...489`
## • `` -> `...490`
## • `` -> `...491`
## • `` -> `...492`
## • `` -> `...493`
## • `` -> `...494`
## • `` -> `...495`
## • `` -> `...496`
## • `` -> `...497`
## • `` -> `...498`
## • `` -> `...499`
## • `` -> `...500`
## • `` -> `...501`
## • `` -> `...502`
## • `` -> `...503`
## • `` -> `...504`
## • `` -> `...505`
## • `` -> `...506`
## • `` -> `...507`
## • `` -> `...508`
## • `` -> `...509`
## • `` -> `...510`
## • `` -> `...511`
## • `` -> `...512`
## • `` -> `...513`
## • `` -> `...514`
## • `` -> `...515`
## • `` -> `...516`
## • `` -> `...517`
## • `` -> `...518`
## • `` -> `...519`
## • `` -> `...520`
## • `` -> `...521`
## • `` -> `...522`
## • `` -> `...523`
## • `` -> `...524`
## • `` -> `...525`
## • `` -> `...526`
## • `` -> `...527`
## • `` -> `...528`
## • `` -> `...529`
## • `` -> `...530`
## • `` -> `...531`
## • `` -> `...532`
## • `` -> `...533`
## • `` -> `...534`
## • `` -> `...535`
## • `` -> `...536`
## • `` -> `...537`
## • `` -> `...538`
## • `` -> `...539`
## • `` -> `...540`
## • `` -> `...541`
## • `` -> `...542`
## • `` -> `...543`
## • `` -> `...544`
## • `` -> `...545`
## • `` -> `...546`
## • `` -> `...547`
## • `` -> `...548`
## • `` -> `...549`
## • `` -> `...550`
## • `` -> `...551`
## • `` -> `...552`
## • `` -> `...553`
## • `` -> `...554`
## • `` -> `...555`
## • `` -> `...556`
## • `` -> `...557`
## • `` -> `...558`
## • `` -> `...559`
## • `` -> `...560`
## • `` -> `...561`
## • `` -> `...562`
## • `` -> `...563`
## • `` -> `...564`
## • `` -> `...565`
## • `` -> `...566`
## • `` -> `...567`
## • `` -> `...568`
## • `` -> `...569`
## • `` -> `...570`
## • `` -> `...571`
## • `` -> `...572`
## • `` -> `...573`
## • `` -> `...574`
## • `` -> `...575`
## • `` -> `...576`
## • `` -> `...577`
## • `` -> `...578`
## • `` -> `...579`
## • `` -> `...580`
## • `` -> `...581`
## • `` -> `...582`
## • `` -> `...583`
## • `` -> `...584`
## • `` -> `...585`
## • `` -> `...586`
## • `` -> `...587`
## • `` -> `...588`
## • `` -> `...589`
## • `` -> `...590`
## • `` -> `...591`
## • `` -> `...592`
## • `` -> `...593`
## • `` -> `...594`
## • `` -> `...595`
## • `` -> `...596`
## • `` -> `...597`
## • `` -> `...598`
## • `` -> `...599`
## • `` -> `...600`
## • `` -> `...601`
## • `` -> `...602`
## • `` -> `...603`
## • `` -> `...604`
## • `` -> `...605`
## • `` -> `...606`
## • `` -> `...607`
## • `` -> `...608`
## • `` -> `...609`
## • `` -> `...610`
## • `` -> `...611`
## • `` -> `...612`
## • `` -> `...613`
## • `` -> `...614`
## • `` -> `...615`
## • `` -> `...616`
## • `` -> `...617`
## • `` -> `...618`
## • `` -> `...619`
## • `` -> `...620`
## • `` -> `...621`
## • `` -> `...622`
## • `` -> `...623`
## • `` -> `...624`
## • `` -> `...625`
## • `` -> `...626`
## • `` -> `...627`
## • `` -> `...628`
## • `` -> `...629`
## • `` -> `...630`
## • `` -> `...631`
## • `` -> `...632`
## • `` -> `...633`
## • `` -> `...634`
## • `` -> `...635`
## • `` -> `...636`
## • `` -> `...637`
## • `` -> `...638`
## • `` -> `...639`
## • `` -> `...640`
## • `` -> `...641`
## • `` -> `...642`
## • `` -> `...643`
## • `` -> `...644`
## • `` -> `...645`
## • `` -> `...646`
## • `` -> `...647`
## • `` -> `...648`
## • `` -> `...649`
## • `` -> `...650`
## • `` -> `...651`
## • `` -> `...652`
## • `` -> `...653`
## • `` -> `...654`
## • `` -> `...655`
## • `` -> `...656`
## • `` -> `...657`
## • `` -> `...658`
## • `` -> `...659`
## • `` -> `...660`
## • `` -> `...661`
## • `` -> `...662`
## • `` -> `...663`
## • `` -> `...664`
## • `` -> `...665`
## • `` -> `...666`
## • `` -> `...667`
## • `` -> `...668`
## • `` -> `...669`
## • `` -> `...670`
## • `` -> `...671`
## • `` -> `...672`
## • `` -> `...673`
## • `` -> `...674`
## • `` -> `...675`
## • `` -> `...676`
## • `` -> `...677`
## • `` -> `...678`
## • `` -> `...679`
## • `` -> `...680`
## • `` -> `...681`
## • `` -> `...682`
## • `` -> `...683`
## • `` -> `...684`
## • `` -> `...685`
## • `` -> `...686`
## • `` -> `...687`
## • `` -> `...688`
## • `` -> `...689`
## • `` -> `...690`
## • `` -> `...691`
## • `` -> `...692`
## • `` -> `...693`
## • `` -> `...694`
## • `` -> `...695`
## • `` -> `...696`
## • `` -> `...697`
## • `` -> `...698`
## • `` -> `...699`
## • `` -> `...700`
## • `` -> `...701`
## • `` -> `...702`
## • `` -> `...703`
## • `` -> `...704`
## • `` -> `...705`
## • `` -> `...706`
## • `` -> `...707`
## • `` -> `...708`
## • `` -> `...709`
## • `` -> `...710`
## • `` -> `...711`
## • `` -> `...712`
## • `` -> `...713`
## • `` -> `...714`
## • `` -> `...715`
## • `` -> `...716`
## • `` -> `...717`
## • `` -> `...718`
## • `` -> `...719`
## • `` -> `...720`
## • `` -> `...721`
## • `` -> `...722`
## • `` -> `...723`
## • `` -> `...724`
## • `` -> `...725`
## • `` -> `...726`
## • `` -> `...727`
## • `` -> `...728`
## • `` -> `...729`
## • `` -> `...730`
## • `` -> `...731`
## • `` -> `...732`
## • `` -> `...733`
## • `` -> `...734`
## • `` -> `...735`
## • `` -> `...736`
## • `` -> `...737`
## • `` -> `...738`
## • `` -> `...739`
## • `` -> `...740`
## • `` -> `...741`
## • `` -> `...742`
## • `` -> `...743`
## • `` -> `...744`
## • `` -> `...745`
## • `` -> `...746`
## • `` -> `...747`
## • `` -> `...748`
## • `` -> `...749`
## • `` -> `...750`
## • `` -> `...751`
## • `` -> `...752`
## • `` -> `...753`
## • `` -> `...754`
## • `` -> `...755`
## • `` -> `...756`
## • `` -> `...757`
## • `` -> `...758`
## • `` -> `...759`
## • `` -> `...760`
## • `` -> `...761`
## • `` -> `...762`
## • `` -> `...763`
## • `` -> `...764`
## • `` -> `...765`
## • `` -> `...766`
## • `` -> `...767`
## • `` -> `...768`
## • `` -> `...769`
## • `` -> `...770`
## • `` -> `...771`
## • `` -> `...772`
## • `` -> `...773`
## • `` -> `...774`
## • `` -> `...775`
## • `` -> `...776`
## • `` -> `...777`
## • `` -> `...778`
## • `` -> `...779`
## • `` -> `...780`
## • `` -> `...781`
## • `` -> `...782`
## • `` -> `...783`
## • `` -> `...784`
## • `` -> `...785`
## • `` -> `...786`
## • `` -> `...787`
## • `` -> `...788`
## • `` -> `...789`
## • `` -> `...790`
## • `` -> `...791`
## • `` -> `...792`
## • `` -> `...793`
## • `` -> `...794`
## • `` -> `...795`
## • `` -> `...796`
## • `` -> `...797`
## • `` -> `...798`
## • `` -> `...799`
## • `` -> `...800`
## • `` -> `...801`
## • `` -> `...802`
## • `` -> `...803`
## • `` -> `...804`
## • `` -> `...805`
## • `` -> `...806`
## • `` -> `...807`
## • `` -> `...808`
## • `` -> `...809`
## • `` -> `...810`
## • `` -> `...811`
## • `` -> `...812`
## • `` -> `...813`
## • `` -> `...814`
## • `` -> `...815`
## • `` -> `...816`
## • `` -> `...817`
## • `` -> `...818`
## • `` -> `...819`
## • `` -> `...820`
## • `` -> `...821`
## • `` -> `...822`
## • `` -> `...823`
## • `` -> `...824`
## • `` -> `...825`
## • `` -> `...826`
## • `` -> `...827`
## • `` -> `...828`
## • `` -> `...829`
## • `` -> `...830`
## • `` -> `...831`
## • `` -> `...832`
## • `` -> `...833`
## • `` -> `...834`
## • `` -> `...835`
## • `` -> `...836`
## • `` -> `...837`
## • `` -> `...838`
## • `` -> `...839`
## • `` -> `...840`
## • `` -> `...841`
## • `` -> `...842`
## • `` -> `...843`
## • `` -> `...844`
## • `` -> `...845`
## • `` -> `...846`
## • `` -> `...847`
## • `` -> `...848`
## • `` -> `...849`
## • `` -> `...850`
## • `` -> `...851`
## • `` -> `...852`
## • `` -> `...853`
## • `` -> `...854`
## • `` -> `...855`
## • `` -> `...856`
## • `` -> `...857`
## • `` -> `...858`
## • `` -> `...859`
## • `` -> `...860`
## • `` -> `...861`
## • `` -> `...862`
## • `` -> `...863`
## • `` -> `...864`
## • `` -> `...865`
## • `` -> `...866`
## • `` -> `...867`
## • `` -> `...868`
## • `` -> `...869`
## • `` -> `...870`
## • `` -> `...871`
## • `` -> `...872`
## • `` -> `...873`
## • `` -> `...874`
## • `` -> `...875`
## • `` -> `...876`
## • `` -> `...877`
## • `` -> `...878`
## • `` -> `...879`
## • `` -> `...880`
## • `` -> `...881`
## • `` -> `...882`
## • `` -> `...883`
## • `` -> `...884`
## • `` -> `...885`
## • `` -> `...886`
## • `` -> `...887`
## • `` -> `...888`
## • `` -> `...889`
## • `` -> `...890`
## • `` -> `...891`
## • `` -> `...892`
## • `` -> `...893`
## • `` -> `...894`
## • `` -> `...895`
## • `` -> `...896`
## • `` -> `...897`
## • `` -> `...898`
## • `` -> `...899`
## • `` -> `...900`
## • `` -> `...901`
## • `` -> `...902`
## • `` -> `...903`
## • `` -> `...904`
## • `` -> `...905`
## • `` -> `...906`
## • `` -> `...907`
## • `` -> `...908`
## • `` -> `...909`
## • `` -> `...910`
## • `` -> `...911`
## • `` -> `...912`
## • `` -> `...913`
## • `` -> `...914`
## • `` -> `...915`
## • `` -> `...916`
## • `` -> `...917`
## • `` -> `...918`
## • `` -> `...919`
## • `` -> `...920`
## • `` -> `...921`
## • `` -> `...922`
## • `` -> `...923`
## • `` -> `...924`
## • `` -> `...925`
## • `` -> `...926`
## • `` -> `...927`
## • `` -> `...928`
## • `` -> `...929`
## • `` -> `...930`
## • `` -> `...931`
## • `` -> `...932`
## • `` -> `...933`
## • `` -> `...934`
## • `` -> `...935`
## • `` -> `...936`
## • `` -> `...937`
## • `` -> `...938`
## • `` -> `...939`
## • `` -> `...940`
## • `` -> `...941`
## • `` -> `...942`
## • `` -> `...943`
## • `` -> `...944`
## • `` -> `...945`
## • `` -> `...946`
## • `` -> `...947`
## • `` -> `...948`
## • `` -> `...949`
## • `` -> `...950`
## • `` -> `...951`
## • `` -> `...952`
## • `` -> `...953`
## • `` -> `...954`
## • `` -> `...955`
## • `` -> `...956`
## • `` -> `...957`
## • `` -> `...958`
## • `` -> `...959`
## • `` -> `...960`
## • `` -> `...961`
## • `` -> `...962`
## • `` -> `...963`
## • `` -> `...964`
## • `` -> `...965`
## • `` -> `...966`
## • `` -> `...967`
## • `` -> `...968`
## • `` -> `...969`
## • `` -> `...970`
## • `` -> `...971`
## • `` -> `...972`
## • `` -> `...973`
## • `` -> `...974`
## • `` -> `...975`
## • `` -> `...976`
## • `` -> `...977`
## • `` -> `...978`
## • `` -> `...979`
## • `` -> `...980`
## • `` -> `...981`
## • `` -> `...982`
## • `` -> `...983`
## • `` -> `...984`
## • `` -> `...985`
## • `` -> `...986`
## • `` -> `...987`
## • `` -> `...988`
## • `` -> `...989`
## • `` -> `...990`
## • `` -> `...991`
## • `` -> `...992`
## • `` -> `...993`
## • `` -> `...994`
## • `` -> `...995`
## • `` -> `...996`
## • `` -> `...997`
## • `` -> `...998`
## • `` -> `...999`
## • `` -> `...1000`
## • `` -> `...1001`
## • `` -> `...1002`
## • `` -> `...1003`
## • `` -> `...1004`
## • `` -> `...1005`
## • `` -> `...1006`
## • `` -> `...1007`
## • `` -> `...1008`
## • `` -> `...1009`
## • `` -> `...1010`
## • `` -> `...1011`
## • `` -> `...1012`
## • `` -> `...1013`
## • `` -> `...1014`
## • `` -> `...1015`
## • `` -> `...1016`
## • `` -> `...1017`
## • `` -> `...1018`
## • `` -> `...1019`
## • `` -> `...1020`
## • `` -> `...1021`
## • `` -> `...1022`
## • `` -> `...1023`
## • `` -> `...1024`
## • `` -> `...1025`
## • `` -> `...1026`
## • `` -> `...1027`
## • `` -> `...1028`
## • `` -> `...1029`
## • `` -> `...1030`
## • `` -> `...1031`
## • `` -> `...1032`
## • `` -> `...1033`
## • `` -> `...1034`
## • `` -> `...1035`
## • `` -> `...1036`
## • `` -> `...1037`
## • `` -> `...1038`
## • `` -> `...1039`
## • `` -> `...1040`
## • `` -> `...1041`
## • `` -> `...1042`
## • `` -> `...1043`
## • `` -> `...1044`
## • `` -> `...1045`
## • `` -> `...1046`
## • `` -> `...1047`
## • `` -> `...1048`
## • `` -> `...1049`
## • `` -> `...1050`
## • `` -> `...1051`
## • `` -> `...1052`
## • `` -> `...1053`
## • `` -> `...1054`
## • `` -> `...1055`
## • `` -> `...1056`
## • `` -> `...1057`
## • `` -> `...1058`
## • `` -> `...1059`
## • `` -> `...1060`
## • `` -> `...1061`
## • `` -> `...1062`
## • `` -> `...1063`
## • `` -> `...1064`
## • `` -> `...1065`
## • `` -> `...1066`
## • `` -> `...1067`
## • `` -> `...1068`
## • `` -> `...1069`
## • `` -> `...1070`
## • `` -> `...1071`
## • `` -> `...1072`
## • `` -> `...1073`
## • `` -> `...1074`
## • `` -> `...1075`
## • `` -> `...1076`
## • `` -> `...1077`
## • `` -> `...1078`
## • `` -> `...1079`
## • `` -> `...1080`
## • `` -> `...1081`
## • `` -> `...1082`
## • `` -> `...1083`
## • `` -> `...1084`
## • `` -> `...1085`
## • `` -> `...1086`
## • `` -> `...1087`
## • `` -> `...1088`
## • `` -> `...1089`
## • `` -> `...1090`
## • `` -> `...1091`
## • `` -> `...1092`
## • `` -> `...1093`
## • `` -> `...1094`
## • `` -> `...1095`
## • `` -> `...1096`
## • `` -> `...1097`
## • `` -> `...1098`
## • `` -> `...1099`
## • `` -> `...1100`
## • `` -> `...1101`
## • `` -> `...1102`
## • `` -> `...1103`
## • `` -> `...1104`
## • `` -> `...1105`
## • `` -> `...1106`
## • `` -> `...1107`
## • `` -> `...1108`
## • `` -> `...1109`
## • `` -> `...1110`
## • `` -> `...1111`
## • `` -> `...1112`
## • `` -> `...1113`
## • `` -> `...1114`
## • `` -> `...1115`
## • `` -> `...1116`
## • `` -> `...1117`
## • `` -> `...1118`
## • `` -> `...1119`
## • `` -> `...1120`
## • `` -> `...1121`
## • `` -> `...1122`
## • `` -> `...1123`
## • `` -> `...1124`
## • `` -> `...1125`
## • `` -> `...1126`
## • `` -> `...1127`
## • `` -> `...1128`
## • `` -> `...1129`
## • `` -> `...1130`
## • `` -> `...1131`
## • `` -> `...1132`
## • `` -> `...1133`
## • `` -> `...1134`
## • `` -> `...1135`
## • `` -> `...1136`
## • `` -> `...1137`
## • `` -> `...1138`
## • `` -> `...1139`
## • `` -> `...1140`
## • `` -> `...1141`
## • `` -> `...1142`
## • `` -> `...1143`
## • `` -> `...1144`
## • `` -> `...1145`
## • `` -> `...1146`
## • `` -> `...1147`
## • `` -> `...1148`
## • `` -> `...1149`
## • `` -> `...1150`
## • `` -> `...1151`
## • `` -> `...1152`
## • `` -> `...1153`
## • `` -> `...1154`
## • `` -> `...1155`
## • `` -> `...1156`
## • `` -> `...1157`
## • `` -> `...1158`
## • `` -> `...1159`
## • `` -> `...1160`
## • `` -> `...1161`
## • `` -> `...1162`
## • `` -> `...1163`
## • `` -> `...1164`
## • `` -> `...1165`
## • `` -> `...1166`
## • `` -> `...1167`
## • `` -> `...1168`
## • `` -> `...1169`
## • `` -> `...1170`
## • `` -> `...1171`
## • `` -> `...1172`
## • `` -> `...1173`
## • `` -> `...1174`
## • `` -> `...1175`
## • `` -> `...1176`
## • `` -> `...1177`
## • `` -> `...1178`
## • `` -> `...1179`
## • `` -> `...1180`
## • `` -> `...1181`
## • `` -> `...1182`
## • `` -> `...1183`
## • `` -> `...1184`
## • `` -> `...1185`
## • `` -> `...1186`
## • `` -> `...1187`
## • `` -> `...1188`
## • `` -> `...1189`
## • `` -> `...1190`
## • `` -> `...1191`
## • `` -> `...1192`
## • `` -> `...1193`
## • `` -> `...1194`
## • `` -> `...1195`
## • `` -> `...1196`
## • `` -> `...1197`
## • `` -> `...1198`
## • `` -> `...1199`
## • `` -> `...1200`
## • `` -> `...1201`
## • `` -> `...1202`
## • `` -> `...1203`
## • `` -> `...1204`
## • `` -> `...1205`
## • `` -> `...1206`
## • `` -> `...1207`
## • `` -> `...1208`
## • `` -> `...1209`
## • `` -> `...1210`
## • `` -> `...1211`
## • `` -> `...1212`
## • `` -> `...1213`
## • `` -> `...1214`
## • `` -> `...1215`
## • `` -> `...1216`
## • `` -> `...1217`
## • `` -> `...1218`
## • `` -> `...1219`
## • `` -> `...1220`
## • `` -> `...1221`
## • `` -> `...1222`
## • `` -> `...1223`
## • `` -> `...1224`
## • `` -> `...1225`
## • `` -> `...1226`
## • `` -> `...1227`
## • `` -> `...1228`
## • `` -> `...1229`
## • `` -> `...1230`
## • `` -> `...1231`
## • `` -> `...1232`
## • `` -> `...1233`
## • `` -> `...1234`
## • `` -> `...1235`
## • `` -> `...1236`
## • `` -> `...1237`
## • `` -> `...1238`
## • `` -> `...1239`
## • `` -> `...1240`
## • `` -> `...1241`
## • `` -> `...1242`
## • `` -> `...1243`
## • `` -> `...1244`
## • `` -> `...1245`
## • `` -> `...1246`
## • `` -> `...1247`
## • `` -> `...1248`
## • `` -> `...1249`
## • `` -> `...1250`
## • `` -> `...1251`
## • `` -> `...1252`
## • `` -> `...1253`
## • `` -> `...1254`
## • `` -> `...1255`
## • `` -> `...1256`
## • `` -> `...1257`
## • `` -> `...1258`
## • `` -> `...1259`
## • `` -> `...1260`
## • `` -> `...1261`
## • `` -> `...1262`
## • `` -> `...1263`
## • `` -> `...1264`
## • `` -> `...1265`
## • `` -> `...1266`
## • `` -> `...1267`
## • `` -> `...1268`
## • `` -> `...1269`
## • `` -> `...1270`
## • `` -> `...1271`
## • `` -> `...1272`
## • `` -> `...1273`
## • `` -> `...1274`
## • `` -> `...1275`
## • `` -> `...1276`
## • `` -> `...1277`
## • `` -> `...1278`
## • `` -> `...1279`
## • `` -> `...1280`
## • `` -> `...1281`
## • `` -> `...1282`
## • `` -> `...1283`
## • `` -> `...1284`
## • `` -> `...1285`
## • `` -> `...1286`
## • `` -> `...1287`
## • `` -> `...1288`
## • `` -> `...1289`
## • `` -> `...1290`
## • `` -> `...1291`
## • `` -> `...1292`
## • `` -> `...1293`
## • `` -> `...1294`
## • `` -> `...1295`
## • `` -> `...1296`
## • `` -> `...1297`
## • `` -> `...1298`
## • `` -> `...1299`
## • `` -> `...1300`
## • `` -> `...1301`
## • `` -> `...1302`
## • `` -> `...1303`
## • `` -> `...1304`
## • `` -> `...1305`
## • `` -> `...1306`
## • `` -> `...1307`
## • `` -> `...1308`
## • `` -> `...1309`
## • `` -> `...1310`
## • `` -> `...1311`
## • `` -> `...1312`
## • `` -> `...1313`
## • `` -> `...1314`
## • `` -> `...1315`
## • `` -> `...1316`
## • `` -> `...1317`
## • `` -> `...1318`
## • `` -> `...1319`
## • `` -> `...1320`
## • `` -> `...1321`
## • `` -> `...1322`
## • `` -> `...1323`
## • `` -> `...1324`
## • `` -> `...1325`
## • `` -> `...1326`
## • `` -> `...1327`
## • `` -> `...1328`
## • `` -> `...1329`
## • `` -> `...1330`
## • `` -> `...1331`
## • `` -> `...1332`
## • `` -> `...1333`
## • `` -> `...1334`
## • `` -> `...1335`
## • `` -> `...1336`
## • `` -> `...1337`
## • `` -> `...1338`
## • `` -> `...1339`
## • `` -> `...1340`
## • `` -> `...1341`
## • `` -> `...1342`
## • `` -> `...1343`
## • `` -> `...1344`
## • `` -> `...1345`
## • `` -> `...1346`
## • `` -> `...1347`
## • `` -> `...1348`
## • `` -> `...1349`
## • `` -> `...1350`
## • `` -> `...1351`
## • `` -> `...1352`
## • `` -> `...1353`
## • `` -> `...1354`
## • `` -> `...1355`
## • `` -> `...1356`
## • `` -> `...1357`
## • `` -> `...1358`
## • `` -> `...1359`
## • `` -> `...1360`
## • `` -> `...1361`
## • `` -> `...1362`
## • `` -> `...1363`
## • `` -> `...1364`
## • `` -> `...1365`
## • `` -> `...1366`
## • `` -> `...1367`
## • `` -> `...1368`
## • `` -> `...1369`
## • `` -> `...1370`
## • `` -> `...1371`
## • `` -> `...1372`
## • `` -> `...1373`
## • `` -> `...1374`
## • `` -> `...1375`
## • `` -> `...1376`
## • `` -> `...1377`
## • `` -> `...1378`
## • `` -> `...1379`
## • `` -> `...1380`
## • `` -> `...1381`
## • `` -> `...1382`
## • `` -> `...1383`
## • `` -> `...1384`
## • `` -> `...1385`
## • `` -> `...1386`
## • `` -> `...1387`
## • `` -> `...1388`
## • `` -> `...1389`
## • `` -> `...1390`
## • `` -> `...1391`
## • `` -> `...1392`
## • `` -> `...1393`
## • `` -> `...1394`
## • `` -> `...1395`
## • `` -> `...1396`
## • `` -> `...1397`
## • `` -> `...1398`
## • `` -> `...1399`
## • `` -> `...1400`
## • `` -> `...1401`
## • `` -> `...1402`
## • `` -> `...1403`
## • `` -> `...1404`
## • `` -> `...1405`
## • `` -> `...1406`
## • `` -> `...1407`
## • `` -> `...1408`
## • `` -> `...1409`
## • `` -> `...1410`
## • `` -> `...1411`
## • `` -> `...1412`
## • `` -> `...1413`
## • `` -> `...1414`
## • `` -> `...1415`
## • `` -> `...1416`
## • `` -> `...1417`
## • `` -> `...1418`
## • `` -> `...1419`
## • `` -> `...1420`
## • `` -> `...1421`
## • `` -> `...1422`
## • `` -> `...1423`
## • `` -> `...1424`
## • `` -> `...1425`
## • `` -> `...1426`
## • `` -> `...1427`
## • `` -> `...1428`
## • `` -> `...1429`
## • `` -> `...1430`
## • `` -> `...1431`
## • `` -> `...1432`
## • `` -> `...1433`
## • `` -> `...1434`
## • `` -> `...1435`
## • `` -> `...1436`
## • `` -> `...1437`
## • `` -> `...1438`
## • `` -> `...1439`
## • `` -> `...1440`
## • `` -> `...1441`
## • `` -> `...1442`
## • `` -> `...1443`
## • `` -> `...1444`
## • `` -> `...1445`
## • `` -> `...1446`
## • `` -> `...1447`
## • `` -> `...1448`
## • `` -> `...1449`
## • `` -> `...1450`
## • `` -> `...1451`
## • `` -> `...1452`
## • `` -> `...1453`
## • `` -> `...1454`
## • `` -> `...1455`
## • `` -> `...1456`
## • `` -> `...1457`
## • `` -> `...1458`
## • `` -> `...1459`
## • `` -> `...1460`
## • `` -> `...1461`
## • `` -> `...1462`
## • `` -> `...1463`
## • `` -> `...1464`
## • `` -> `...1465`
## • `` -> `...1466`
## • `` -> `...1467`
## • `` -> `...1468`
## • `` -> `...1469`
## • `` -> `...1470`
## • `` -> `...1471`
## • `` -> `...1472`
## • `` -> `...1473`
## • `` -> `...1474`
## • `` -> `...1475`
## • `` -> `...1476`
## • `` -> `...1477`
## • `` -> `...1478`
## • `` -> `...1479`
## • `` -> `...1480`
## • `` -> `...1481`
## • `` -> `...1482`
## • `` -> `...1483`
## • `` -> `...1484`
## • `` -> `...1485`
## • `` -> `...1486`
## • `` -> `...1487`
## • `` -> `...1488`
## • `` -> `...1489`
## • `` -> `...1490`
## • `` -> `...1491`
## • `` -> `...1492`
## • `` -> `...1493`
## • `` -> `...1494`
## • `` -> `...1495`
## • `` -> `...1496`
## • `` -> `...1497`
## • `` -> `...1498`
## • `` -> `...1499`
## • `` -> `...1500`
## • `` -> `...1501`
## • `` -> `...1502`
## • `` -> `...1503`
## • `` -> `...1504`
## • `` -> `...1505`
## • `` -> `...1506`
## • `` -> `...1507`
## • `` -> `...1508`
## • `` -> `...1509`
## • `` -> `...1510`
## • `` -> `...1511`
## • `` -> `...1512`
## • `` -> `...1513`
## • `` -> `...1514`
## • `` -> `...1515`
## • `` -> `...1516`
## • `` -> `...1517`
## • `` -> `...1518`
## • `` -> `...1519`
## • `` -> `...1520`
## • `` -> `...1521`
## • `` -> `...1522`
## • `` -> `...1523`
## • `` -> `...1524`
## • `` -> `...1525`
## • `` -> `...1526`
## • `` -> `...1527`
## • `` -> `...1528`
## • `` -> `...1529`
## • `` -> `...1530`
## • `` -> `...1531`
## • `` -> `...1532`
## • `` -> `...1533`
## • `` -> `...1534`
## • `` -> `...1535`
## • `` -> `...1536`
## • `` -> `...1537`
## • `` -> `...1538`
## • `` -> `...1539`
## • `` -> `...1540`
## • `` -> `...1541`
## • `` -> `...1542`
## • `` -> `...1543`
## • `` -> `...1544`
## • `` -> `...1545`
## • `` -> `...1546`
## • `` -> `...1547`
## • `` -> `...1548`
## • `` -> `...1549`
## • `` -> `...1550`
## • `` -> `...1551`
## • `` -> `...1552`
## • `` -> `...1553`
## • `` -> `...1554`
## • `` -> `...1555`
## • `` -> `...1556`
## • `` -> `...1557`
## • `` -> `...1558`
## • `` -> `...1559`
## • `` -> `...1560`
## • `` -> `...1561`
## • `` -> `...1562`
## • `` -> `...1563`
## • `` -> `...1564`
## • `` -> `...1565`
## • `` -> `...1566`
## • `` -> `...1567`
## • `` -> `...1568`
## • `` -> `...1569`
## • `` -> `...1570`
## • `` -> `...1571`
## • `` -> `...1572`
## • `` -> `...1573`
## • `` -> `...1574`
## • `` -> `...1575`
## • `` -> `...1576`
## • `` -> `...1577`
## • `` -> `...1578`
## • `` -> `...1579`
## • `` -> `...1580`
## • `` -> `...1581`
## • `` -> `...1582`
## • `` -> `...1583`
## • `` -> `...1584`
## • `` -> `...1585`
## • `` -> `...1586`
## • `` -> `...1587`
## • `` -> `...1588`
## • `` -> `...1589`
## • `` -> `...1590`
## • `` -> `...1591`
## • `` -> `...1592`
## • `` -> `...1593`
## • `` -> `...1594`
## • `` -> `...1595`
## • `` -> `...1596`
## • `` -> `...1597`
## • `` -> `...1598`
## • `` -> `...1599`
## • `` -> `...1600`
## • `` -> `...1601`
## • `` -> `...1602`
## • `` -> `...1603`
## • `` -> `...1604`
## • `` -> `...1605`
## • `` -> `...1606`
## • `` -> `...1607`
## • `` -> `...1608`
## • `` -> `...1609`
## • `` -> `...1610`
## • `` -> `...1611`
## • `` -> `...1612`
## • `` -> `...1613`
## • `` -> `...1614`
## • `` -> `...1615`
## • `` -> `...1616`
## • `` -> `...1617`
## • `` -> `...1618`
## • `` -> `...1619`
## • `` -> `...1620`
## • `` -> `...1621`
## • `` -> `...1622`
## • `` -> `...1623`
## • `` -> `...1624`
## • `` -> `...1625`
## • `` -> `...1626`
## • `` -> `...1627`
## • `` -> `...1628`
## • `` -> `...1629`
## • `` -> `...1630`
## • `` -> `...1631`
## • `` -> `...1632`
## • `` -> `...1633`
## • `` -> `...1634`
## • `` -> `...1635`
## • `` -> `...1636`
## • `` -> `...1637`
## • `` -> `...1638`
## • `` -> `...1639`
## • `` -> `...1640`
## • `` -> `...1641`
## • `` -> `...1642`
## • `` -> `...1643`
## • `` -> `...1644`
## • `` -> `...1645`
## • `` -> `...1646`
## • `` -> `...1647`
## • `` -> `...1648`
## • `` -> `...1649`
## • `` -> `...1650`
## • `` -> `...1651`
## • `` -> `...1652`
## • `` -> `...1653`
## • `` -> `...1654`
## • `` -> `...1655`
## • `` -> `...1656`
## • `` -> `...1657`
## • `` -> `...1658`
## • `` -> `...1659`
## • `` -> `...1660`
## • `` -> `...1661`
## • `` -> `...1662`
## • `` -> `...1663`
## • `` -> `...1664`
## • `` -> `...1665`
## • `` -> `...1666`
## • `` -> `...1667`
## • `` -> `...1668`
## • `` -> `...1669`
## • `` -> `...1670`
## • `` -> `...1671`
## • `` -> `...1672`
## • `` -> `...1673`
## • `` -> `...1674`
## • `` -> `...1675`
## • `` -> `...1676`
## • `` -> `...1677`
## • `` -> `...1678`
## • `` -> `...1679`
## • `` -> `...1680`
## • `` -> `...1681`
## • `` -> `...1682`
## • `` -> `...1683`
## • `` -> `...1684`
## • `` -> `...1685`
## • `` -> `...1686`
## • `` -> `...1687`
## • `` -> `...1688`
## • `` -> `...1689`
## • `` -> `...1690`
## • `` -> `...1691`
## • `` -> `...1692`
## • `` -> `...1693`
## • `` -> `...1694`
## • `` -> `...1695`
## • `` -> `...1696`
## • `` -> `...1697`
## • `` -> `...1698`
## • `` -> `...1699`
## • `` -> `...1700`
## • `` -> `...1701`
## • `` -> `...1702`
## • `` -> `...1703`
## • `` -> `...1704`
## • `` -> `...1705`
## • `` -> `...1706`
## • `` -> `...1707`
## • `` -> `...1708`
## • `` -> `...1709`
## • `` -> `...1710`
## • `` -> `...1711`
## • `` -> `...1712`
## • `` -> `...1713`
## • `` -> `...1714`
## • `` -> `...1715`
## • `` -> `...1716`
## • `` -> `...1717`
## • `` -> `...1718`
## • `` -> `...1719`
## • `` -> `...1720`
## • `` -> `...1721`
## • `` -> `...1722`
## • `` -> `...1723`
## • `` -> `...1724`
## • `` -> `...1725`
## • `` -> `...1726`
## • `` -> `...1727`
## • `` -> `...1728`
## • `` -> `...1729`
## • `` -> `...1730`
## • `` -> `...1731`
## • `` -> `...1732`
## • `` -> `...1733`
## • `` -> `...1734`
## • `` -> `...1735`
## • `` -> `...1736`
## • `` -> `...1737`
## • `` -> `...1738`
## • `` -> `...1739`
## • `` -> `...1740`
## • `` -> `...1741`
## • `` -> `...1742`
## • `` -> `...1743`
## • `` -> `...1744`
## • `` -> `...1745`
## • `` -> `...1746`
## • `` -> `...1747`
## • `` -> `...1748`
## • `` -> `...1749`
## • `` -> `...1750`
## • `` -> `...1751`
## • `` -> `...1752`
## • `` -> `...1753`
## • `` -> `...1754`
## • `` -> `...1755`
## • `` -> `...1756`
## • `` -> `...1757`
## • `` -> `...1758`
## • `` -> `...1759`
## • `` -> `...1760`
## • `` -> `...1761`
## • `` -> `...1762`
## • `` -> `...1763`
## • `` -> `...1764`
## • `` -> `...1765`
## • `` -> `...1766`
## • `` -> `...1767`
## • `` -> `...1768`
## • `` -> `...1769`
## • `` -> `...1770`
## • `` -> `...1771`
## • `` -> `...1772`
## • `` -> `...1773`
## • `` -> `...1774`
## • `` -> `...1775`
## • `` -> `...1776`
## • `` -> `...1777`
## • `` -> `...1778`
## • `` -> `...1779`
## • `` -> `...1780`
## • `` -> `...1781`
## • `` -> `...1782`
## • `` -> `...1783`
## • `` -> `...1784`
## • `` -> `...1785`
## • `` -> `...1786`
## • `` -> `...1787`
## • `` -> `...1788`
## • `` -> `...1789`
## • `` -> `...1790`
## • `` -> `...1791`
## • `` -> `...1792`
## • `` -> `...1793`
## • `` -> `...1794`
## • `` -> `...1795`
## • `` -> `...1796`
## • `` -> `...1797`
## • `` -> `...1798`
## • `` -> `...1799`
## • `` -> `...1800`
## • `` -> `...1801`
## • `` -> `...1802`
## • `` -> `...1803`
## • `` -> `...1804`
## • `` -> `...1805`
## • `` -> `...1806`
## • `` -> `...1807`
## • `` -> `...1808`
## • `` -> `...1809`
## • `` -> `...1810`
## • `` -> `...1811`
## • `` -> `...1812`
## • `` -> `...1813`
## • `` -> `...1814`
## • `` -> `...1815`
## • `` -> `...1816`
## • `` -> `...1817`
## • `` -> `...1818`
## • `` -> `...1819`
## • `` -> `...1820`
## • `` -> `...1821`
## • `` -> `...1822`
## • `` -> `...1823`
## • `` -> `...1824`
## • `` -> `...1825`
## • `` -> `...1826`
## • `` -> `...1827`
## • `` -> `...1828`
## • `` -> `...1829`
## • `` -> `...1830`
## • `` -> `...1831`
## • `` -> `...1832`
## • `` -> `...1833`
## • `` -> `...1834`
## • `` -> `...1835`
## • `` -> `...1836`
## • `` -> `...1837`
## • `` -> `...1838`
## • `` -> `...1839`
## • `` -> `...1840`
## • `` -> `...1841`
## • `` -> `...1842`
## • `` -> `...1843`
## • `` -> `...1844`
## • `` -> `...1845`
## • `` -> `...1846`
## • `` -> `...1847`
## • `` -> `...1848`
## • `` -> `...1849`
## • `` -> `...1850`
## • `` -> `...1851`
## • `` -> `...1852`
## • `` -> `...1853`
## • `` -> `...1854`
## • `` -> `...1855`
## • `` -> `...1856`
## • `` -> `...1857`
## • `` -> `...1858`
## • `` -> `...1859`
## • `` -> `...1860`
## • `` -> `...1861`
## • `` -> `...1862`
## • `` -> `...1863`
## • `` -> `...1864`
## • `` -> `...1865`
## • `` -> `...1866`
## • `` -> `...1867`
## • `` -> `...1868`
## • `` -> `...1869`
## • `` -> `...1870`
## • `` -> `...1871`
## • `` -> `...1872`
## • `` -> `...1873`
## • `` -> `...1874`
## • `` -> `...1875`
## • `` -> `...1876`
## • `` -> `...1877`
## • `` -> `...1878`
## • `` -> `...1879`
## • `` -> `...1880`
## • `` -> `...1881`
## • `` -> `...1882`
## • `` -> `...1883`
## • `` -> `...1884`
## • `` -> `...1885`
## • `` -> `...1886`
## • `` -> `...1887`
## • `` -> `...1888`
## • `` -> `...1889`
## • `` -> `...1890`
## • `` -> `...1891`
## • `` -> `...1892`
## • `` -> `...1893`
## • `` -> `...1894`
## • `` -> `...1895`
## • `` -> `...1896`
## • `` -> `...1897`
## • `` -> `...1898`
## • `` -> `...1899`
## • `` -> `...1900`
## • `` -> `...1901`
## • `` -> `...1902`
## • `` -> `...1903`
## • `` -> `...1904`
## • `` -> `...1905`
## • `` -> `...1906`
## • `` -> `...1907`
## • `` -> `...1908`
## • `` -> `...1909`
## • `` -> `...1910`
## • `` -> `...1911`
## • `` -> `...1912`
## • `` -> `...1913`
## • `` -> `...1914`
## • `` -> `...1915`
## • `` -> `...1916`
## • `` -> `...1917`
## • `` -> `...1918`
## • `` -> `...1919`
## • `` -> `...1920`
## • `` -> `...1921`
## • `` -> `...1922`
## • `` -> `...1923`
## • `` -> `...1924`
## • `` -> `...1925`
## • `` -> `...1926`
## • `` -> `...1927`
## • `` -> `...1928`
## • `` -> `...1929`
## • `` -> `...1930`
## • `` -> `...1931`
## • `` -> `...1932`
## • `` -> `...1933`
## • `` -> `...1934`
## • `` -> `...1935`
## • `` -> `...1936`
## • `` -> `...1937`
## • `` -> `...1938`
## • `` -> `...1939`
## • `` -> `...1940`
## • `` -> `...1941`
## • `` -> `...1942`
## • `` -> `...1943`
## • `` -> `...1944`
## • `` -> `...1945`
## • `` -> `...1946`
## • `` -> `...1947`
## • `` -> `...1948`
## • `` -> `...1949`
## • `` -> `...1950`
## • `` -> `...1951`
## • `` -> `...1952`
## • `` -> `...1953`
## • `` -> `...1954`
## • `` -> `...1955`
## • `` -> `...1956`
## • `` -> `...1957`
## • `` -> `...1958`
## • `` -> `...1959`
## • `` -> `...1960`
## • `` -> `...1961`
## • `` -> `...1962`
## • `` -> `...1963`
## • `` -> `...1964`
## • `` -> `...1965`
## • `` -> `...1966`
## • `` -> `...1967`
## • `` -> `...1968`
## • `` -> `...1969`
## • `` -> `...1970`
## • `` -> `...1971`
## • `` -> `...1972`
## • `` -> `...1973`
## • `` -> `...1974`
## • `` -> `...1975`
## • `` -> `...1976`
## • `` -> `...1977`
## • `` -> `...1978`
## • `` -> `...1979`
## • `` -> `...1980`
## • `` -> `...1981`
## • `` -> `...1982`
## • `` -> `...1983`
## • `` -> `...1984`
## • `` -> `...1985`
## • `` -> `...1986`
## • `` -> `...1987`
## • `` -> `...1988`
## • `` -> `...1989`
## • `` -> `...1990`
## • `` -> `...1991`
## • `` -> `...1992`
## • `` -> `...1993`
## • `` -> `...1994`
## • `` -> `...1995`
## • `` -> `...1996`
## • `` -> `...1997`
## • `` -> `...1998`
## • `` -> `...1999`
## • `` -> `...2000`
## • `` -> `...2001`
## • `` -> `...2002`
## • `` -> `...2003`
## • `` -> `...2004`
## • `` -> `...2005`
## • `` -> `...2006`
## • `` -> `...2007`
## • `` -> `...2008`
## • `` -> `...2009`
## • `` -> `...2010`
## • `` -> `...2011`
## • `` -> `...2012`
## • `` -> `...2013`
## • `` -> `...2014`
## • `` -> `...2015`
## • `` -> `...2016`
## • `` -> `...2017`
## • `` -> `...2018`
## • `` -> `...2019`
## • `` -> `...2020`
## • `` -> `...2021`
## • `` -> `...2022`
## • `` -> `...2023`
## • `` -> `...2024`
## • `` -> `...2025`
## • `` -> `...2026`
## • `` -> `...2027`
## • `` -> `...2028`
## • `` -> `...2029`
## • `` -> `...2030`
## • `` -> `...2031`
## • `` -> `...2032`
## • `` -> `...2033`
## • `` -> `...2034`
## • `` -> `...2035`
## • `` -> `...2036`
## • `` -> `...2037`
## • `` -> `...2038`
## • `` -> `...2039`
## • `` -> `...2040`
## • `` -> `...2041`
## • `` -> `...2042`
## • `` -> `...2043`
## • `` -> `...2044`
## • `` -> `...2045`
## • `` -> `...2046`
## • `` -> `...2047`
## • `` -> `...2048`
## • `` -> `...2049`
## • `` -> `...2050`
## • `` -> `...2051`
## • `` -> `...2052`
## • `` -> `...2053`
## • `` -> `...2054`
## • `` -> `...2055`
## • `` -> `...2056`
## • `` -> `...2057`
## • `` -> `...2058`
## • `` -> `...2059`
## • `` -> `...2060`
## • `` -> `...2061`
## • `` -> `...2062`
## • `` -> `...2063`
## • `` -> `...2064`
## • `` -> `...2065`
## • `` -> `...2066`
## • `` -> `...2067`
## • `` -> `...2068`
## • `` -> `...2069`
## • `` -> `...2070`
## • `` -> `...2071`
## • `` -> `...2072`
## • `` -> `...2073`
## • `` -> `...2074`
## • `` -> `...2075`
## • `` -> `...2076`
## • `` -> `...2077`
## • `` -> `...2078`
## • `` -> `...2079`
## • `` -> `...2080`
## • `` -> `...2081`
## • `` -> `...2082`
## • `` -> `...2083`
## • `` -> `...2084`
## • `` -> `...2085`
## • `` -> `...2086`
## • `` -> `...2087`
## • `` -> `...2088`
## • `` -> `...2089`
## • `` -> `...2090`
## • `` -> `...2091`
## • `` -> `...2092`
## • `` -> `...2093`
## • `` -> `...2094`
## • `` -> `...2095`
## • `` -> `...2096`
## • `` -> `...2097`
## • `` -> `...2098`
## • `` -> `...2099`
## • `` -> `...2100`
## • `` -> `...2101`
## • `` -> `...2102`
## • `` -> `...2103`
## • `` -> `...2104`
## • `` -> `...2105`
## • `` -> `...2106`
## • `` -> `...2107`
## • `` -> `...2108`
## • `` -> `...2109`
## • `` -> `...2110`
## • `` -> `...2111`
## • `` -> `...2112`
## • `` -> `...2113`
## • `` -> `...2114`
## • `` -> `...2115`
## • `` -> `...2116`
## • `` -> `...2117`
## • `` -> `...2118`
## • `` -> `...2119`
## • `` -> `...2120`
## • `` -> `...2121`
## • `` -> `...2122`
## • `` -> `...2123`
## • `` -> `...2124`
## • `` -> `...2125`
## • `` -> `...2126`
## • `` -> `...2127`
## • `` -> `...2128`
## • `` -> `...2129`
## • `` -> `...2130`
## • `` -> `...2131`
## • `` -> `...2132`
## • `` -> `...2133`
## • `` -> `...2134`
## • `` -> `...2135`
## • `` -> `...2136`
## • `` -> `...2137`
## • `` -> `...2138`
## • `` -> `...2139`
## • `` -> `...2140`
## • `` -> `...2141`
## • `` -> `...2142`
## • `` -> `...2143`
## • `` -> `...2144`
## • `` -> `...2145`
## • `` -> `...2146`
## • `` -> `...2147`
## • `` -> `...2148`
## • `` -> `...2149`
## • `` -> `...2150`
## • `` -> `...2151`
## • `` -> `...2152`
## • `` -> `...2153`
## • `` -> `...2154`
## • `` -> `...2155`
## • `` -> `...2156`
## • `` -> `...2157`
## • `` -> `...2158`
## • `` -> `...2159`
## • `` -> `...2160`
## • `` -> `...2161`
## • `` -> `...2162`
## • `` -> `...2163`
## • `` -> `...2164`
## • `` -> `...2165`
## • `` -> `...2166`
# define and solve prioritization problem
prob <- problem(units, pfeatures) %>%
      add_max_features_objective(budget = prot_area + 10) %>%
      add_feature_weights(weights = ds$tree$edge.length) %>%
      add_relative_targets(0.1) %>%
      add_binary_decisions() %>%
      add_locked_in_constraints(locked_in = protected) %>%
      add_rsymphony_solver(gap = 0)
soln <- solve(prob)
plot(soln + protected, 
     col = c("gray", "orange", "darkgreen"))

Exercises:

  • These exercises were based on the chronogram, which we packaged up with the California flora dataset. Replicate them with the phylogram and cladogram. Do the conservation priorities change?
  • prioritizr has functionality to consider many additional factors, such as land cost, reserve shape, and reserve connectivity. A simple example is add_neighbor_constraints(k = 2), which would force proposed reserves to have at least two neighbors. Like other constraints this will likely add a lot of computational time, so try it when you have a few minutes to wait.
prob <- problem(units, features) %>%
      add_relative_targets(0.1) %>%
      add_binary_decisions() %>%
      add_rsymphony_solver(gap = 0)

soln1 <- prob %>%
      add_max_phylo_div_objective(budget = prot_area,
                                  tree = chronogram) %>%
      solve()
soln2 <- prob %>%
      add_max_phylo_div_objective(budget = prot_area,
                                  tree = cladogram) %>%
      solve()
soln3 <- prob %>%
      add_max_phylo_div_objective(budget = prot_area,
                                  tree = phylogram) %>%
      solve(run_checks = F)

solns <- stack(soln1, soln2, soln3)

# plot priorities for each facet
solns %>%
      setNames(c("chronogram", "cladogram", "phylogram")) %>%
      plot(main = "cladogram",
           col = c("gray", "red"))

# plot number of facets under which each cell is a priority
solns %>%
      sum() %>%
      plot(col = c("gray", "yellow", "orange", "red"))


And that is all!